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Abstract 

A time accurate, high-order, conservative, yet efficient method named Finite Spectral Volume (FSV) is 
developed for conservation laws on unstructured grids . The concept of a “spectral volume’ is introduced 
to achieve high-order accuracy in an efficient manner similar to spectral element and multi-domain 
spectral methods. In addition, each spectral volume is further sub-divided into control volumes (CVs), 
and cell-averaged data from these control volumes is used to reconstruct a high-order approximation in 
the spectral volume. Riemann solvers are used to compute the fluxes at spectral volume boundaries. Then 
cell-averaged state variables in the control volumes are updated independently. Furthermore, TVD (Total 
Variation Diminishing) and TVB (Total Variation Bounded) limiters are introduced in the FSV method to 
remove/reduce spurious oscillations near discontinuities. A very desirable feature of the F SV method is 
that the reconstruction is carried out only once, and analytically, and is the same for all cells of the same 
type, and that the reconstruction stencil is always non-singular, in contrast to the memory and CPU- 
intensive reconstruction in a high-order finite volume (FV) method. Discussions are made concerning 
why the FSV method is significantly more efficient than high-order finite volume and the Discontinuous 
Galerkin (DG) methods. Fundamental properties of the FSV method are studied and high-order accuracy 
is demonstrated for several model problems with and without discontinuities. 

1 Framework of the Finite Spectral Volume Method 

To present the basic idea, we consider the following multi-dimensional scalar conservation laws: 

du(x<,y y t) df (u{x, y, f)) dg(u(x, y,0) __ ^ 

— 5T~ + F + F 0 (U) 

on domain Q with the following initial condition 

M(x,y,0) = u Q (x y y) (lb) 

and appropriate boundary conditions on. Domain Q is discretized into N non-overlapping cells which are 
called spectral volumes ( SVs ), i.e. 

N 

Q = (JS, (2) 

i=l 

The reason why the cells are called SVs will be clear later. Integrating (1) on a SV S,, we obtain 



0 


where F = (f g), and n is the unit outward normal of 3S, , the boundary of S r Define the cell averaged 
state variable for S, as 

J udV 

U, = ^77— (4) 


where V, is the volume (area in 2D) of S,-. Then (3) becomes 

^- + — \ {(F»n)dA = 0 (5) 

dt V 

i r=l A, 

where L is the number of faces (edges in 2D) in S„ and A r represents the r-th face. The surface integration 
on each face can be performed with a k-th order accurate Gauss quadrature formula, i.e. 


J(F • tl)clA = Yj W rj F ( U ( X rj.y r j )) • n r A r + ^(4^ ) 


where \v r j are the Gauss quadrature weights, (x rp y rj ) are the Gauss quadrature points, h is the maximum 
span of all the SVs in x and y directions, time t is omitted whenever there is no confusion. If F — constant , 
the following identity exists: 

SJ ( F»n)dA = 0 (7) 

r=1 4 

Therefore, we will gain an extra order if we sum up the surface integrals for all faces of S„ i.e., 

X \(F *n)dA= ^^w rj F(u(x r j y rj )) • it r A r + 0(A r h k + l ) (8) 

r=l A r r-\j=\ 

Since 0(V, ) = 0(A r h ) , therefore we have 

* n)dA^^-Y J ^w rj F(u(.x rj y r j))*n r A r +0(h k ) (9) 

v i r=l Ar V <r=\j = \ 

Now assume a multi-dimensional polynomial in x and y of order at most k - 1 exists on 5, which is a k-th 
order approximation to the state variable, i.e., 

p i (x,y) = u(x,y) + 0(h k ), (x,y)eS, (10) 

This polynomial is called a reconstruction of the state variable. With the polynomial distribution on each 
SV , the state variable is most likely discontinuous across the SV boundaries, unless the state variable is a 
polynomial of order k - 1 or less. Therefore, the flux integration involves two discontinuous state 
variables just to the left and right of a face of the SV boundary. This flux integration is then carried out 
using an exact Riemann solver or one of the Lipschitz continuous approximate Riemann solvers or flux 
splitting procedures, i.e., 

F(u{x j > y rj )) 9 ^ Riemann ( Pi * Try )’ Pn (*rj ’ T rj )’ ) 

j j (I1 ) 

+ 0( pi (x rj , y rj ) - p n (x rj , y rj )) 

Here p n is the reconstruction polynomial of a neighboring SV S m which shares face A r with 5,. Both p { 
and p n are k-th order approximations of the exact state variable, i.e., 

Pi (x rj ,y r j) = ll i x rj » y r j ) + °( f,k ) ( 1 2a ) 

Pn ( x ri - >V/ ) = u ( x ri - T rj ) + 0 i hk ) ( 1 2b ) 


Therefore 
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F(u(x r j , y r j )) • n = F Rjemann (pj(x rj , y rj ), p n (x rj , y rj ), n r ) + 0(h k ) (13) 

Substituting (13) into (6), we obtain 

J 

j(F »n)dA- Yj w rj F Riema,m(Pi(x rj ,y rj ), p u (x rj , y rj ),n r )A r + 0(A r h ) ( 14 ) 

\ j =l 

Summarizing (5)-(14). we obtain the following semi-discrete, k-th order accurate scheme on S, for the 
conservation laws (1) 


du, 1 

— y- H 

dt V. 


L J 

11 - 


rj r Riemann 


r=l 7 =l 


( Pi ('^Vy ’ y rj Pn rj ’ 3 rj) ’ 0(h ) 


(15) 


What we have done so far follows exactly the finite volume doctrine. We, however, omitted a vital detail, 
i.e., how we build the high-order reconstruction polynomial given just the cell-averaged state variables for 
the SVs. Here is where the new method departs from the traditional FV scheme. In a FV method, a stencil 
(a group of neighboring cells and the cell under consideration) is used to build a high-order polynomial 
approximation to the state variable on the cell under consideration. Depending on the order of accuracy, a 
very large number of (up to 60-100) cells may be necessary to perform a non-singular quadratic data 
reconstruction. For an unstructured grid, each cell has a unique reconstruction stencil, and the 
reconstruction problem needs to be solved for each and every cell at each and every time step or iteration. 
The reconstruction can be very memory and CPU intensive especially for higher than linear 
reconstructions. This is probably why we have seen few attempts to perform quadratic reconstructions in 
three dimensions. A recent study on a 3-D quadratic reconstruction [5] showed that the cost in memory 
and CPU time does not justify the effort. 


In this paper, the FSV method is developed to address this very drawback. Our solution is as follows. 
Instead of using a large stencil of neighboring cells to perform the reconstruction, we subdivide the SV 
into control volumes (CVV). The order of accuracy of the reconstruction determines the number of CVs to 
be generated in each SV. For example, for a linear reconstruction on a triangle, the triangle is divided into 
at least three CVs as shown in Figure la, and cell averaged state variables are defined on the CVs . Figures 
lb and lc give some possible CV subdivisions for quadratic and cubic data reconstructions. The number 
of CVs in Figure 1 is the minimum required for these polynomial reconstructions. Other CV subdivisions 
are definitely possible. 


With any of these high-order reconstructions, Eq. (15) can then be used to update the cell-averaged state 
variable on the SVs , i.e., the cell-averaged state variable w, for S t at a new time level n+7 (i.e. u * 1 ) can 
be obtained with an appropriate time integration scheme based on the solution at time level n with t = 
nAt , where At is the time step. However, in order to use the same high-order reconstruction at time level 

n+1 , it is necessary to “scatter” the update A u i = — u" back to the cell-averaged state variables at 

all the CVs in S,. This is how we perform the scattering operation. Each CV inside a SV is treated 
separately os if it is independent to update the cell-averaged state variable for the CV. Note that the subtle 
difference between a FV and a FSV method is that all the CVs in a SV use the same data reconstruction. 
As a result, it is not necessary to use a Riemann flux or flux splitting for the interior boundaries between 
the CVs inside a particular SV because the state variable is continuous across the interior CV boundaries. 
Riemann fluxes are only necessary at the boundaries of the SV. To maintain a high-order accuracy, 
Gaussian quadrature formulas are again used, not only for the Riemann fluxes through the SV boundaries, 
but also for the fluxes through interior CV boundaries. The most significant advantage of the FSV method, 
as compared with the FV method, is that the reconstruction for a particular cell type (e.g. triangles) with a 
certain CV subdivision (e.g. those shown in Figure 1) is exactly the same. Even though the shape of the 
SVs may all be different, as long as they are geometric similarly subdivided, they all have the same 
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reconstruction (in the parameter space, to be explained in the final paper), and the weights tor evaluating 
the state variables in term of cell-averaged unknowns at similar quadrature points are all the same. 
Therefore, the memory and CPU intensive reconstructions used in a FV method are solved analytically 
without taking any extra memory in the FSV method. Furthermore, exact fluxes rather than Riemann 
fluxes are used at the interior boundaries of the CVs, resulting again significant savings because the 
Riemann flux is usually several time more expensive to compute than the exact flux. 

The idea can of course be easily extended to other cell types such as quadrilaterals, tetrahedra, hexahedra, 
prisms, etc. For cell types other than triangles and tetrahedra, it appears that symmetric CV subdivisions 
with the minimum number of CVs for a given order of accuracy are not easily found. 

The FSV method shares many advantages with the Discontinuous Galerkin (DG) method [2-4] in that it is 
compact which is suitable for parallel computing, high-order accurate, conservative, and capable of 
handling complex geometries. Furthermore, the FSV method is expected to be much more efficient than 
the DG method, and has higher resolution than the DG method for discontinuities because of the 
availability of local cell-averaged state variables at the control volumes. 

The main steps in a k-th order FSV method (with an order k-1 polynomial reconstruction) are: 

1. Compute the state variables at the quadrature points; 

2. Use a k-th order accurate quadrature formula (exact for a polynomial of order k-1) and a Riemann 
solver to compute the surface flux integrals at the spectral volume boundaries, and use a k-th order 
accurate quadrature formula for analytical fluxes for interior control volume boundaries because the 
state-variable is continuous across the interior CV boundaries; 

3. Use a TVD Runge-Kutta scheme for time integration; 

The main steps in a k-th order DG method are: 

1. Compute the state variables at the quadrature points; 

2. Use a Riemann solver and a 2k-th order quadrature formula to compute the surface flux integrals; 

3. Use a ( 2k-l)-th order quadrature formula to compute the volume integrals; 

4. Left multiply the residual by the inverse of the mass matrix because the mass matrix is usually not 
diagonal for k > 2; 

5. Use a TVD Runge-Kutta scheme for time integration; 

Note that with the FSV method, the high-order volume integration required in a DG method is completely 
eliminated. Furthermore, the surface integral in the FSV method needs only to be k-th order accurate 
instead of the 2k-th order accuracy required in a DG method. As a result, the FSV method requires only 
half the quadrature points required by a DG method to carry out the surface integration. For fourth order 
DG and FSV methods, there are 10 degrees of freedom (DOF) in 2D for a single variable. In a DG 
method, 4 quadrature points are required to compute the surface integral on a single edge to achieve the 
desired accuracy, and 12 quadrature points are required to compute the volume integral up to the desired 
accuracy [2]. To update all the DOFs (assuming a single variable) for a single element using the DG 
method, (3x4x10 + 12x10) = 240 variable evaluations at the quadrature points are required. In order to 
compute the surface integrals, 3x4x10 =120 Riemann fluxes need to be computed. In contrast, only 2 
quadrature points are necessary to compute the surface integral on a single edge to achieve the desired 
accuracy with a fourth-order FSV method. Therefore, to update all the DOFs for a single spectral volume 
(element) using the FSV method, only 2x27 = 54 (27 being the number of total edges in the spectral 
volume shown in Figure lc) variable evaluations at the quadrature points are required. In addition, only 
2x12 = 24 (12 being the total number of spectral volume boundary edges) Riemann fluxes are required, 
and the rest 2x15 = 30 (15 being the number of interior control volume boundary-edges) fluxes are 
analytical fluxes because the reconstructed state variable is continuous across the interior control volume 
boundaries inside the spectral volume. It is well known that a Riemann flux is usually several times more 
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expensive to compute than the analytical flux for the Euler equations. Let's assume that a Riemann flux is 
only three times as expensive as an analytical flux (a very conservative estimate indeed). Then the FSV 
method requires 24+30/3 = 34 Riemann fluxes. The third-order TVD Runge-Kutta scheme takes 
negligible CPU time because of the very few number of operations. Finally in a fourth DG method, the 
residual vector has to left-multiply by a 10x10 matrix at each iteration. 

If we assume that the Riemann flux computation dominates the total CPU time, then the FSV method is 
about 4 times as fast as the DG method. If on the other hand, variable evaluations dominate the CPU time, 
the FSV method can be close to five times as fast as the DG method. Overall, we expect the fourth-order 
FSV method to be about 4-5 times as fast as the fourth-order DG method in 2D. If one is interested in 
even higher order accurate schemes, the availability of very high-order quadrature formulas may become 
an issue in a DG method. For example, a sixth order DG scheme necessitates a 12 th order quadrature 
formula for surface integration, and an ll Lh order quadrature formula for volume integration. In three- 
dimensions, it is expected the advantage of the FSV method is even more pronounced because high-order 
quadrature formulas for volume integration in a tetrahedron are required in the DG method. 

2. Numerical Tests 

We have implemented the FSV method in both ID and 2D, with a variety of limiters (control-volume- 
wise (CV-wise) and spectral-volume-wise (SV-wise) TVDM and TVBM limiters) to eliminate spurious 
oscillations. In the final paper, detailed formulations will be given. Here we just show several 
demonstration cases to demonstrate its capability. 


1 Test with the Burger's Equation 

In this test, we solve the non-linear Burger’s equation with a periodic boundary condition: 

du du 2 / 2 , 

dt dx 


u(x, 0) = iiq(x) = 1 + — sin(7ix). 

The exact solution is smooth up to t = 2/n , then it develops a moving shock, which interacts with 
rarefaction waves. At t = 0.3, the solution is still smooth. FSV schemes from second to sixth order of 
accuracy are tested, and the Lj and L Ioo errors are listed in Table 1, together with the numerical order of 
accuracy. Note that the expected formal orders of accuracy for all the tested schemes are achieved in both 
the Lj and L <OQ norms. The computed solution with a second-order FSV scheme on 6 SVs is compared 
with the solution with a fourth-order FSV scheme on 3 SVs in Figure 2. The numerical solutions therefore 
have the same number of degrees-of-freedom. Note that the fourth-order scheme gave a visibly better 
solution than the second-order scheme. 


At t = 2/n, a shock starts to form in the solution. The numerical solution would be oscillatory without 
limiters. Figure 3 displays the computed solutions with a 4 th -order FSV scheme on 20 SVs using various 
limiters. Note that the SVTVDM limiter strongly dissipated the numerical solution, while the CVTVDM 
limiter gave a much better solution. Both TVBM limiters with M = 20 gave reasonable results, with the 
CVTVBM limiter showing a slightly more accurate prediction. The solution with the SVTVBM limiter is 
slightly oscillatory. 

At t = 1, a shock wave has formed in the solution. The numerical solutions computed with a fourth-order 
FSV scheme on 20 and 40 SVs using both TVBM limiters are presented in Figure 4. Note that the shock- 
wave is generally captured in one spectral volume, and the CVTVBM limiter once again produced a 
solution with a better resolution for the shock wave. 
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Finally to see whether TVBM limiters affect the solution accuracy away from the shock wave, the local 
solution errors with a 4 th -order FSV scheme on different grids are plotted in Figure 5. It is shown clearly 
in the Figure that the expected order of accuracy (4 th -order) is retained away from the shock-wave with 
both the CVTVBM and SVTVBM limiters. The CVTVBM limiter is shown again to have a better 
resolution for the shock wave. 


2. Accuracy Study with 2D Linear Wave Equation 

In this case, we test the accuracy of the FSV method on the linear equation: 

du du da „ , 

— + — + — = 0, - 1 < jc < 1, - 1 < v < I 

dt dx dy 

u(x , y,0) = Uq(x, y), periodic boundary condition 


The initial condition is iiq(x, y) = sin 7t(x + y ) . The numerical simulation is carried until t = 1 on a 
triangular grid generated from a uniform Cartesian grid by cutting each Cartesian cell into two triangles. 
In Table 2, we show the Lj and L , TO errors produced using second to fourth order FSV schemes with CVs 
shown in Figure 1. The third-order TVD Runge-Kutta time-integration scheme was used for all the 
computations presented here. The errors presented in the table are time-step independent because the time 
step At was made small enough so that the errors are dominated by the spatial discretization. Again it is 
shown that the desired order of accuracy is obtained for all cases. 


Conclusions 

A high-order Finite Spectral Volume (FSV) method is developed in this study. The concept of “spectral 
volume” is introduced to achieve high-order accuracy in a very efficient manner. The FSV method is 
much more efficient in terms of both memory and CPU requirement than a high-order finite volume 
method because the reconstruction for a particular grid type is solved only once, and analytically, and is 
never explicitly carried out. Furthermore, the “reconstruction stencils” are never singular. We also 
discussed why the FSV method is significantly more efficient than the DG method. Control-volume-wise 
and spectral-volume-wise TVDM and TVBM limiters are developed to remove spurious oscillations near 
discontinuities. It has been shown that CV-wise limiters perform better than SV-wise limiters. Because of 
the availability of local data, the FSV method is expected to produce much sharper discontinuity profiles 
than the DG method. 


Accuracy studies with ID and 2D linear and non-linear scalar conservation laws have been carried out, 
and the order of accuracy claim has been numerically verified. The TVBM limiters were found to 
maintain uniformly high-order accuracy away from discontinuities. 
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Table 2. Accuracy on u, +u x +« v = 0, u 0 (x,y) = sin[7r(x + y)] at t = 1 


Order of Accuracy 

Grid 

Ly error 

Lj order 

error 

order 

2 

10x10x2 

1.64e-2 

- 

4.13e-2 

- 

20x20x2 

4.01e-3 

2.03 

9.59e-3 

2.11 

40x40x2 

9.85e-4 

2.03 

2.21e-3 

2.18 

80x80x2 

2.44e-4 

2.01 

5.18e-4 

2.09 

160x160x2 

6.09e-5 

2.00 

1.24e-4 

2.06 

3 

10x10x2 

4.18e-3 

- 

7.76e-3 

- 

20x20x2 

5.33e-4 

2.97 

1.01e-3 

2.94 

40x40x2 

6.70e-5 

2.99 

1.24e-4 

3.03 

80x80x2 

8.13e-6 

3.04 

1 .5 le-5 

3.04 

160x160x2 

1.05e-6 

2.95 

1.93e-6 

2.97 

4 

10x10x2 

9.33e-5 

- 

3.17e-4 

- 

20x20x2 

5.86e-6 

3.99 

1.94e-5 

4.02 

40x40x2 

3.70e-7 

3.99 

1.24e-6 

3.97 

80x80x2 

2.32e-8 

4.00 

7.78e-8 

3.99 

160x160x2 

1.45e-9 

4.00 

4.84e-9 

4.01 
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Figure 1. Control Volumes in a Triangular Spectral Volume (a) Linear; (b) Quadratic; (c) Cubic. 



Figure 2. Computed Solutions to the Burger's Equation at t = 0.3 Using a Second-Order and Fourth Order 
FSV Schemes with 12 Degrees-of-Freedom without Limiters 
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Figure 3. Computed Solutions to the Burger's Equation at t = 2 /ti Using a Fourth Order FSV Scheme with 
CVTVDM, SVTVDM, and CVTVBM, SVTVBM Limiters on 10 SVs 



Figure 4. Computed Solutions to the Burger’s Equations at t = 1 with CVTVBM and SVTVBM Limiters 

Using 20 and 40 Spectral Volumes 



Figure 5. Local Error of Computed Solutions of the Burger's Equation at t = 1 with a Fourth-Order FSV 

scheme and CVTVBM and SVTVBM Limiters 
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